// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                      avtSphericalCompactnessFactorQuery.C                 //
// ************************************************************************* //

#include <avtSphericalCompactnessFactorQuery.h>

#include <vtkCell.h>

#ifdef PARALLEL
#include <mpi.h>
#endif

#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkVisItUtility.h>

#include <avtCallback.h>
#include <avtParallel.h>
#include <avtRevolvedVolume.h>
#include <avtSourceFromAVTDataset.h>
#include <avtVMetricVolume.h>

#include <DebugStream.h>
#include <ImproperUseException.h>

#include <float.h>
#include <math.h>

using     std::string;
using     std::vector;

#if !defined(M_PI)
#define M_PI 3.14159265358979323846
#endif

// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery constructor
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
// ****************************************************************************

avtSphericalCompactnessFactorQuery::avtSphericalCompactnessFactorQuery()
{
    rev_volume = new avtRevolvedVolume;
    volume = new avtVMetricVolume;
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery destructor
//
//  Purpose:
//      Defines the destructor.  Note: this should not be inlined in the header
//      because it causes problems for certain compilers.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
// ****************************************************************************

avtSphericalCompactnessFactorQuery::~avtSphericalCompactnessFactorQuery()
{
    delete volume;
    delete rev_volume;
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::SetInputParams
//
//  Purpose:
//    Set the input parameters.
//
//  Programmer: Cyrus Harrison
//  Creation: Wed Jul 16 15:58:54 PDT 2014
//
//  Modifications:
//
// ****************************************************************************

void
avtSphericalCompactnessFactorQuery::SetInputParams(const MapNode &params)
{
    if(params.HasNumericVectorEntry("centroid"))
    {
        overrideCentroid = true;

        doubleVector cvals;
        params.GetEntry("centroid")->ToDoubleVector(cvals);
        centroid[0] = cvals[0];
        centroid[1] = cvals[1];
        centroid[2] = cvals[2];
    }else
    {
      overrideCentroid = false;
    }
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::PreExecute
//
//  Purpose:
//      This is called before all of the domains are executed.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
//  Modifications:
//    Jeremy Meredith, Thu Feb 15 11:55:03 EST 2007
//    Call inherited PreExecute before everything else.
//
//    Cyrus Harrison, Wed Jul 16 15:52:57 PDT 2014
//    Added support for user selected center.
//
// ****************************************************************************

void
avtSphericalCompactnessFactorQuery::PreExecute(void)
{
    avtTwoPassDatasetQuery::PreExecute();

    if(!overrideCentroid)
    {
        centroid[0] = 0.;
        centroid[1] = 0.;
        centroid[2] = 0.;
    }

    total_volume = 0.;
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::MidExecute
//
//  Purpose:
//      This is called after Execute1 and before Execute2.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
//  Modifications:
//    Cyrus Harrison, Wed Jul 16 15:52:57 PDT 2014
//    Added support for user selected center.
//
// ****************************************************************************

void
avtSphericalCompactnessFactorQuery::MidExecute(void)
{
    SumDoubleAcrossAllProcessors(total_volume);
    if(!overrideCentroid)
    {
        double C_tmp[3] = { 0, 0, 0 };
        SumDoubleArrayAcrossAllProcessors(centroid, C_tmp, 3);
        if (total_volume != 0.)
        {
            C_tmp[0] /= total_volume;
            C_tmp[1] /= total_volume;
            C_tmp[2] /= total_volume;
        }
        centroid[0] = C_tmp[0];
        centroid[1] = C_tmp[1];
        centroid[2] = C_tmp[2];
    }

    volume_inside = 0.;
    radius = pow(total_volume*0.75/M_PI, 0.3333333);
    if (is2D)
    {
        sphere_center[0] = centroid[0];
        sphere_center[1] = 0.;
        sphere_center[2] = 0.;
    }
    else
    {
        sphere_center[0] = centroid[0];
        sphere_center[1] = centroid[1];
        sphere_center[2] = centroid[2];
    }
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::PostExecute
//
//  Purpose:
//      This is called after all of the domains are executed.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
//  Modifications:
//
//    Cyrus Harrison, Tue Sep 18 13:45:35 PDT 2007
//    Added support for user settable floating point format string
//
//    Kathleen Biagas, Tue Feb 25 16:16:13 PST 2014
//    Add Xml results.
//
//    Cyrus Harrison, Wed Jul 16 15:52:57 PDT 2014
//    Added support for user selected center.
//
// ****************************************************************************

void
avtSphericalCompactnessFactorQuery::PostExecute(void)
{
    SumDoubleAcrossAllProcessors(volume_inside);
    if (total_volume == 0.)
    {
        total_volume = 1;
    }

    //
    //  Parent class uses this message to set the Results message
    //  in the Query Attributes that is sent back to the viewer.
    //  That is all that is required of this query.
    //
    char msg[4096];
    string floatFormat = queryAtts.GetFloatFormat();
    string format = "Spherical Compactness Factor = " + floatFormat 
                    + ".  Using centroid for sphere origin.  Centroid used "
                      "was (" + floatFormat + ", "
                              + floatFormat + ", "
                              + floatFormat + ")"
                      "  Radius was " + floatFormat;
    snprintf(msg, 4096,format.c_str(),
                       volume_inside / total_volume, 
                       sphere_center[0], sphere_center[1], sphere_center[2],
                       radius);
    SetResultMessage(msg);
    SetResultValue(volume_inside / total_volume);
    MapNode result_node;
    result_node["radius"] = radius;

    doubleVector sc(3);
    sc[0] = sphere_center[0];
    sc[1] = sphere_center[1];
    sc[2] = sphere_center[2];

    result_node["override_centroid"] = overrideCentroid;
    result_node["centroid"] = sc;
    result_node["spherical_compactness_factor"] = volume_inside / total_volume;
    SetXmlResult(result_node.ToXML());
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::Execute1
//
//  Purpose:
//      Processes a single domain.  This will calculate the total volume and
//      the centroid on this pass.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
//  Modifications:
//    Kathleen Biagas, Tue Apr 22 07:54:11 MST 2014
//    Use double instead of float.
//
// ****************************************************************************

void
avtSphericalCompactnessFactorQuery::Execute1(vtkDataSet *ds, const int dom)
{
    int nCells = ds->GetNumberOfCells();
    vtkDataArray *ghosts = ds->GetCellData()->GetArray("avtGhostZones");
    vtkDataArray *var = ds->GetCellData()->GetArray("avt_volume");
    if (var == NULL)
    {
        EXCEPTION0(ImproperUseException);
    }
    for (int i = 0 ; i < nCells ; i++)
    {
        if (ghosts != NULL && ghosts->GetTuple1(i) != 0.)
            continue;
        vtkCell *cell = ds->GetCell(i);
        double center[3];
        vtkVisItUtility::GetCellCenter(cell, center);
        double volume = var->GetTuple1(i);
        volume = (volume < 0 ? -volume : volume);

        if(!overrideCentroid)
        {
            centroid[0] += volume*center[0];
            centroid[1] += volume*center[1];
            centroid[2] += volume*center[2];
        }

        total_volume += volume;
    }
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::Execute2
//
//  Purpose:
//      Processes a single domain.  This will calculate the portion inside
//      the sphere on this pass.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
//  Modifications:
//    Kathleen Biagas, Tue Apr 22 07:54:11 MST 2014
//    Use double instead of float.
//
// ****************************************************************************

void
avtSphericalCompactnessFactorQuery::Execute2(vtkDataSet *ds, const int dom)
{
    int nCells = ds->GetNumberOfCells();
    vtkDataArray *ghosts = ds->GetCellData()->GetArray("avtGhostZones");
    vtkDataArray *var = ds->GetCellData()->GetArray("avt_volume");
    if (var == NULL)
    {
        EXCEPTION0(ImproperUseException);
    }
    double rad_squared = radius*radius;
    for (int i = 0 ; i < nCells ; i++)
    {
        if (ghosts != NULL && ghosts->GetTuple1(i) != 0.)
            continue;
        vtkCell *cell = ds->GetCell(i);
        double center[3];
        vtkVisItUtility::GetCellCenter(cell, center);
        double dist = (center[0]-sphere_center[0])*(center[0]-sphere_center[0])
                   + (center[1]-sphere_center[1])*(center[1]-sphere_center[1])
                   + (center[2]-sphere_center[2])*(center[2]-sphere_center[2]);
        if (dist > rad_squared)
            continue;  // Not inside.
        double volume = var->GetTuple1(i);
        volume = (volume < 0 ? -volume : volume);
        volume_inside += volume;
    }
}


// ****************************************************************************
//  Method: avtSphericalCompactnessFactorQuery::ApplyFilters
//
//  Purpose:
//      Applies the filters to the input.
//
//  Programmer: Hank Childs
//  Creation:   July 14, 2005
//
//  Modifications:
//    Kathleen Bonnell, Wed Apr  2 10:44:18 PDT 2008
//    Removed unsued section of code that retrieved varname, retrieved dims
//    and set 'useVar'.
//
// ****************************************************************************

avtDataObject_p
avtSphericalCompactnessFactorQuery::ApplyFilters(avtDataObject_p inData)
{
    //
    // Create an artificial pipeline.
    //
    avtDataset_p ds;
    CopyTo(ds, inData);
    avtSourceFromAVTDataset termsrc(ds);
    avtDataObject_p dob = termsrc.GetOutput();

    //
    // Set up our base class so it is ready to sum.
    //

    avtExpressionFilter *vf = NULL;
    is2D = (dob->GetInfo().GetAttributes().GetSpatialDimension() != 3);
    if (is2D)
        vf = rev_volume;
    else
        vf = volume;

    vf->SetOutputVariableName("avt_volume");
    vf->SetInput(dob);
    dob = vf->GetOutput();
 
    //
    // Cause our artificial pipeline to execute.
    //
    avtContract_p contract =
        inData->GetOriginatingSource()->GetGeneralContract();
    dob->Update(contract);
    return dob;
}


